Conceptual Model
Contents
# import some utilities
# definition of some plotting tools
%matplotlib inline
from bokeh.layouts import gridplot
from bokeh.layouts import layout
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import ColumnDataSource
from bokeh.plotting import output_file
from bokeh.palettes import Category10
from bokeh.io import export_png
from bokeh.models import Legend
from bokeh.models import Range1d
from bokeh.models import ColorBar
from bokeh.models import LinearAxis
from bokeh.models import DatetimeTickFormatter
# select a palette
#from bokeh.palettes import Dark2_25 as palette
from bokeh.palettes import Spectral11, Category20b
from bokeh.palettes import Spectral6
from bokeh.plotting import output_file
from bokeh.transform import linear_cmap
import itertools
import datetime
import itertools
import numpy
import os
import re
import sys
from netCDF4 import Dataset
import pandas
import logging
import math
def color_gen():
yield from itertools.cycle(Category10[10])
color = color_gen()
from bokeh.resources import INLINE
output_notebook(INLINE)
export = False
import logging
import sys
sys.path.append(os.path.join(os.environ["HOME"],"Projects","SeaFlect","python"))
from DataPoolAccessModule import AeolusDataClass
from DataPoolAccessModule import QualityControlClass
from SurfaceReflectanceModule import SurfaceReflectanceModelClass
def make_plot(title, hist, edges):
p = figure(title=title, tools='', background_fill_color="#fafafa")
p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
fill_color="navy", line_color="white", alpha=0.5)
p.y_range.start = 0
#p.legend.location = "center_right"
#p.legend.background_fill_color = "#fefefe"
p.xaxis.axis_label = 'x'
p.yaxis.axis_label = 'Pr(x)'
p.grid.grid_line_color="white"
return p
verbose_level = 'info'
logger=logging.getLogger()
logging.basicConfig(level=getattr(logging, verbose_level.upper()))
Cm = {}
Li = {}
wind_speed = numpy.linspace(2, 24, 23)
theory=SurfaceReflectanceModelClass(logger=logger)
theory.Ro_fixed=0.08
Li["low"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()
Cm["low"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
theory.Ro_fixed=0.016
Li["high"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()
Cm["high"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
AeolusProductObject = AeolusDataClass(debug=False, \
generate=False,
logger=logger)
exp = "fixed_range"
#exp = "flexible_range"
AeolusProductObject.pickle_data_file=os.path.join("SeaFlectProjectData","DataPool","seaflect_data_"+exp + ".pkl")
figure_path = os.path.join("project_figures", "notebook", exp)
if not os.path.exists(figure_path):
os.makedirs(figure_path)
[data_pool, available_periods, available_regions] = AeolusProductObject.get_data()
AeolusProductObject.data_pool=data_pool
# now that we have a consistent dataset we can make an annual
AeolusProductObject._create_annual_(periods=available_periods)
# and globa
AeolusProductObject._create_global_(regions=available_regions)
#definitions of range bins accoring to c-convention
range_bin={"top":0, "atmosphere":3, "surface":4, "ocean":5, "background":6}
Conceptual Model¶
To understand the radiative properties of the composite layer, a simple conceptual model has been developed. This conceptual model considers the macro physical properties of the layer, and link the incoming radiation at the upper and lower interfaces of the layer to the the emerging radiation at these interfaces.
This notebook space here is too short to provide all of the equations, but in words the conceptual model indicates that the upwelling radiation at the top, is the sum of the upwelling radiation at the bottom of the layer multiplied by the transmission of the layer, and the downwelling radiation at the top of the layer multiplied by the reflection of the layer, were we ignore any internal source of radiation. As an example the upwelling radiation at the upper interface of a layer:
where the upper script gives the direction + for upwelling and - for down welling , and the sub-script the location 0 the lower interface, and 1 the upper interface.
A similar equation can be written for the emerging radiation at the lower interface.
The equations are prepared with the assumption that these macro radiative properties \(R\) and \(T\) can be linked to the micro physical properties of a homogeneous layer (e.g. Rayleigh scattering).
With the notion that the composite range bin is composed of a homogeneous atmospheric, surface and oceanograpical layer, we are interested in the macro radiative properties of the composite layer. This can be derived on the basis of the above equations, assuming that there is conservation of energy at the interfaces. The net result is that the macro physical reflection (indicated as \(R_{aso}\), were the a stands for atmosphere, s for surface and o for ocean) can be expressed in terms of these properties for the individual layers.
The second term on right hand side is the multiple reflection term. What we see from this is that the reflection of the composite layer is the reflection of the combined surface and ocean layer (which is what Li et al in their paper considers), extended with a multiple scattering term, as the radiation emerging from the surface-ocean composite layer needs to travel through the atmosphere.
In the end we are interested in an estimation of the \(R_{so}\) term, which can only be derived if we have an estimation of \(R_a\) and \(T_a\).
We can estimate \(R_a\) through the notion that in the clear atmosphere the return signal originates from Rayleigh scattering, which is a function of the number of molucules. We can estimate the Rayleigh scattering from the signal in the range bin above the surface range bin, by taking the number of molecules, into account. Which means we need to normalise signal 1 by the layer thickness in pressure space and multiply this by the pressure difference between the surface pressure and pressure level of the first interface. As we do not have the pressure values, we take the altitudes.
Thus:
were \(S\) represent the measured signal, \(\Delta z\) is the geometrical thickness, and \(z\) is the altitude.
Here we assume a cloud and aerosol free atmospheric layer. It is thus assumed that the transmission of the atmospheric layer is determined by the Rayleigh scattering only. We have no real estimation for this yet (we are building this model), we can fix it to a specific value. (0.9 in the calculations presented here)
In the end we are interested in the \(R_{so}\) as a function of the observatables. Inversion of the conceptual model indicate
it should be noted that \( \left( R_{aso} -R_{a} \right) \) is the estimation of \(R_{so}\) used by Li et al (2010). The second term (\(\left( T^2_a + R_a\left[R_{aso} - R_{a}\right] \right)\)) is a multiple scattering term which accounts for the transmission through the atmospheric layer.
Unfortunately we would need access to the incident radiation for each layer as the reflectance can be estimated as
and we measure the return signal of the layer \(F^+_1\), hence if the incident radiation is known then we can calculate the reflectance, using the auxilliary information to estimate the transmission would close the equations.
The incident radiation at the top of the surface layer is not known, and hence a proxy value for is selected. The uv_total_energy with a value of 60 is a factor 1000 too small as experiments shown. Hence a proxy value of 60000 is used in some cases, but this has no value other than bringing results within a reasonable scale.
region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Raso"][:], color="black", \
legend_label="Composite")
p.yaxis.axis_label = 'Raso [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
t=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="ACCD {} {}".format(period, region))
t.circle(time, dataset["Product"]["Ra"][:], color="red", legend_label="Atmosphere")
t.yaxis.axis_label = 'Ra [-]'
t.y_range=Range1d(-0.05, 0.2)
t.legend.location = "top_left"
graphic.append(t)
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Cm"][:], color="green", legend_label="CM")
p.yaxis.axis_label = 'Reflectance [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Li"][:], color="blue", legend_label="Li")
p.yaxis.axis_label = 'Reflectance [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
figure_name = os.path.join(figure_path,"reflectance_signal_{}_{}.png".format(period, region))
export_png(plot, filename=figure_name)
region="E"
period="AnnualCycle"
qc_level = "basic"
#
AeolusProductObject.snr_up=20
AeolusProductObject.snr_low=10
AeolusProductObject.sunalt_threshold=-7
qc = AeolusProductObject.apply_quality_control(level=qc_level, period=period, region=region)
dataset =AeolusProductObject.data_pool[period][region]
if qc == None:
time_range=slice(kwargs["time_range"][0], kwargs["time_range"][1])
else:
time_range=qc[0]
time = dataset["L1Bmie"]["date"][time_range]
x = dataset["MET"]["surface_wind_speed"][time_range]
#basic signal
col = dataset["L1Bmie"]["snr"][time_range,AeolusProductObject.layer_index["surface"]]
mapper = linear_cmap(field_name='col', palette=Spectral11 ,low=min(col) ,high=max(col))
# cm results
y = dataset["Product"]["Reflectance"]["Cm"][time_range]
source = ColumnDataSource({'x':x,'y':y,'col':col})
graphic=[]
p=figure(plot_width=400, plot_height=300, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="Raw input {} {}".format(period, region))
p.circle( "x", "y", source=source ,line_color=mapper,color=mapper,size = 5)
p.line(wind_speed, Cm["high"][0,:],color="blue", legend_label="high" )
p.line(wind_speed, Cm["low"][0,:],color="green", legend_label="low" )
p.xaxis.axis_label = 'wind speed [m/s]'
p.yaxis.axis_label = 'ACCD '
p.x_range = Range1d(0, 25)
p.y_range=Range1d(0,0.05)
bar = ColorBar(color_mapper=mapper['transform'], width=8, location=(0,0))
p.add_layout(bar, "left")
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=500, toolbar_location=None)
show(plot)
if export:
figure_name = os.path.join(figure_path,"correlation_accd_windspeed_{}_{}.png".format(period, region))
export_png(plot, filename=figure_name)
Summary¶
Based on the analysis presented above, it appears there is a hint of a relation between the wind speed and the surface reflection (first order effect is there). The analysis is performed for the ACCD counts from the Mie channel. The open issues at this point of the analysis are
Conversion of ACCD counts to physical space. This is normally refered to as calibration and we need guidance on this. Our plans are to perform a ROM radiometric calibration using the observations over other radiometrically stable targets.
QC, there is an issue with the implementation of QC, there are many more points passing QC in 2019 than in 2020. This requires further investigation. Also we need to adjust some thresholds to consolidate the results.
Mie channel: should this analysis also be performed for the Rayleigh channel.
A conversion from the Corrected Counts to wind speed according to a Geophysical Model Function (GMF) can only be performed after we have a calibrated signal. We expect however, that the GMF has some weak point with respect to the description of the wind induced white cap fraction, and the reflectance of the white caps. To improve this would require a more detailed analysis using state of the art RT models as the adopted method is based on analysis of photographical material. Such an detailed analysis is beyond scope of current project.
Specific questions¶
how is background measured
is the acc counts corrected for background
how is scattering ratio defined
what if we apply to Rayleigh
would a vicarious calibration effort make sense (there is an ESA ITT out, but not sure it fits there)
calibration equation: what is meaning of Kmie? Values are very large for detector efficiency etc.
what could be reason for degradation of return signal over time (before and after calibration)
We plan to redo this for the data at the higher spatial resolution. Replace atmospheric noise by instrument noise. The accumulation to the product resolution is currently done without taking any properties (atmosphere or dead pixels) into account?
We plan to include dark ocean surface. no results yet.
Do you have a RT code which can be shared and used for a vicarious calibration? we propose to use SmartMC as we anticipate a need for a MC approach. Plan would be to first calibrate the background signal and then port the experience to lidar return.